REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments 

regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 

Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA,  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  oenalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

New  Reprint 


4.  TITLE  AND  SUBTITLE 

Gluing  bifurcations  in  coupled  spin  torque  nano-oscillators 


3.  DATES  COVERED  (From  -  To) 


5a.  CONTRACT  NUMBER 
W91  INF-12-1-0283 


5b.  GRANT  NUMBER 


6.  AUTHORS 

Katherine  Beauvais,  Richard  Shaffer,  Antonio  Palacios,  Visarath  In, 
Teresa  Emery,  Patrick  Longhini,  James  Turtle 


5c.  PROGRAM  ELEMENT  NUMBER 
611102 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

San  Diego  State  University  Foundation 
5250  Campanile  Dr. 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


San  Diego,  CA 


92182  -1931 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND 
ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.Box  12211 

Research  Triangle  Park,  NC  27709-2211 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

60768-EG.2 


12.  DISTRIBUTION  AVAILIBILITY  STATEMENT 
Approved  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


14.  ABSTRACT 

Over  the  past  few  years  it  has  been  shown,  through  theory  and  experiments,  that  the  AC  current  produced  by  spin 
torque  nano-oscillators  (STNO),  coupled  in  an  array,  can  lead  to  feedback  between  the  STNOs  causing  them  to 
synchronize  and  that,  collectively,  the  microwave  power  output  of  the  array  is  significantly  larger  than  that  of  an 
individual  valve.  Other  works  have  pointed,  however,  to  the  difficulty  in  achieving  synchronization.  The  governing 
equations  lead  to  a  highl -dimensional  dynamical  system  with  S_N  symmetry,  i.e.,  permutation  symmetry  of  N 


15.  SUBJECT  TERMS 
Nano-oscillators,  symmetry,  bifurcations 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

15.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

ABSTRACT 

OF  PAGES 

Antonio  Palacios 

UU 

UU 

UU 

UU 

19b.  TELEPHONE  NUMBER 

619-594-6808 

Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39. 1 8 


Report  Title 

Gluing  bifurcations  in  coupled  spin  torque  nano-oscillators 

ABSTRACT 

Over  the  past  few  years  it  has  been  shown,  through  theory  and  experiments,  that  the  AC  current  produced  by  spin 
torque  nano-oscillators  (STNO),  coupled  in  an  array,  can  lead  to  feedback  between  the  STNOs  causing  them  to 
synchronize  and  that,  collectively,  the  microwave  power  output  of  the  array  is  significantly  larger  than  that  of  an 
individual  valve.  Other  works  have  pointed,  however,  to  the  difficulty  in  achieving  synchronization.  The  governing 
equations  lead  to  a  highl-dimensional  dynamical  system  with  S_N  symmetry,  i.e.,  permutation  symmetry  of  N  objects 
due  to  all-to-all  coupling.  Under  the  ARO  grant,  we  have  investigated  a  reduction  of  the  modeling  equations,  via 
center  Manifold  and  Normal  form,  and  the  possibility  of  applying  the  most  recent  results  from  related  work  on 
coupled  Hamiltonian  systems  to  explain  the  nature  of  the  bifurcations. 


REPORT  DOCUMENTATION  PAGE  (SF298) 
(Continuation  Sheet) 


Continuation  for  Block  13 


ARO  Report  Number  60768. 2-EG 
Gluing  bifurcations  in  coupled  spin  torque  nano- 


Block  13:  Supplementary  Note 

©2013  .  Published  in  Journal  of  Applied  Physics,  Vol.  Ed.  0  113,  (11)  (2013),  ( (11).  DoD  Components  reserve  a  royalty-free, 
nonexclusive  and  irrevocable  right  to  reproduce,  publish,  or  otherwise  use  the  work  for  Federal  purposes,  and  to  authroize 
others  to  do  so  (DODGARS  §32.36).  The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and 
should  not  be  construed  as  an  official  Department  of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other 
documentation. 


Approved  for  public  release;  distribution  is  unlimited. 


JOURNAL  OF  APPLIED  PHYSICS  113,  1 14901  (2013) 


®CrossMark 

4-dick  for  updates 


Gluing  bifurcations  in  coupled  spin  torque  nano-oscillators 

James  Turtle,1  Katherine  Beauvais,1  Richard  Shaffer,1  Antonio  Palacios, 1,a)  Visarath  In,2 
Teresa  Emery,2  and  Patrick  Longhini2 

^Nonlinear  Dynamical  Systems  Group,  Department  of  Mathematics,  San  Diego  State  University, 

San  Diego,  California  92182,  USA. 

2Space  and  Naval  Warfare  Systems  Center  Pacific,  Code  71730,  53560  Hull  Street,  San  Diego, 

California  92152-5001,  USA. 

(Received  10  October  2012;  accepted  27  February  2013;  published  online  15  March  2013) 

Over  the  past  few  years,  it  has  been  shown,  through  theory  and  experiments,  that  the  AC  current 
produced  by  spin  torque  nano-oscillators  (STNO),  coupled  in  an  array,  can  lead  to  feedback 
between  the  STNOs  causing  them  to  synchronize  and  that,  collectively,  the  microwave  power 
output  of  the  array  is  significantly  larger  than  that  of  an  individual  valve.  Other  works  have 
pointed,  however,  to  the  difficulty  in  achieving  synchronization.  In  particular,  Persson  et  al. 
[J.  Appl.  Phys.  101,  09A503  (2007)]  shows  that  the  region  of  parameter  space  where  the 
synchronization  state  exists  for  even  a  small  array  with  two  STNOs  is  rather  small.  In  this  work, 
we  explore  in  more  detail  the  nature  of  the  bifurcations  that  lead  into  and  out  of  the 
synchronization  state  for  the  two-array  case.  The  bifurcation  analysis  shows  bistability  between 
in-phase  and  out-of-phase  limit  cycle  oscillations.  In  fact,  there  are  two  distinct  pairs  of  such 
cycles.  But  as  the  input  current  increases,  the  limit  cycles  may  increase  their  amplitudes  until  they 
merge  with  one  another  in  a  gluing  bifurcation.  More  importantly,  we  show  that  changing  the 
direction  of  the  applied  magnetic  held  can,  in  principle,  increase  the  region  of  synchronized 
oscillations.  ©2013  American  Institute  of  Physics.  [http://dx.doi.org/10.1063/L4795266] 


I.  INTRODUCTION 

Spintronics — the  emerging  science  that  seeks  to  exploit 
the  intrinsic  spin  of  the  electron — has  stimulated  scientists 
and  engineers  around  the  world  to  envision,  design,  and  fab¬ 
ricate  an  entire  new  generation  of  smaller,  faster,  and  more 
energy-efficient  nano-electronic  devices.1  Spintronic  devices 
work  on  the  quantum  mechanical  effects  of  electrons  having 
two-state  spins,  “up”  or  “down.”  By  running  current  through 
a  ferromagnetic  material,  a  spin-polarized  current  can  be  cre¬ 
ated  and  manipulated  by  magnetic  fields.  The  most  common 
application  of  this  effect  is  the  spin  nano-valve  device,  which 
consists  of  al  least  two  layers  (about  lOOnm  in  lateral  size) 
of  ferromagnetic  materials  separated  by  a  nonmagnetic  mate¬ 
rial  layer,  see  Fig.  1 .  In  one  layer,  the  magnetization  vectors 
are  fixed  while  on  the  other  hand  they  are  free  in  order  to 
exploit  the  giant  magnetoresistive  (GMR)  effect.  This  effect 
is  observed  as  a  significant  change  in  the  electrical  resistance 
of  some  materials  depending  on  whether  the  magnetization 
of  adjacent  ferromagnetic  layers  is  in  antiparallel  (high  re¬ 
sistance)  or  in  parallel  (low  resistance).  One  immediate 
application  of  the  spin  valve  is  as  a  sensor  of  weak  fields. 
But  a  later  discovery  of  the  spin-polarized  phenomenon  may 
soon  allow  spin  valves  to  be  used  also  as  miniaturized  micro- 
wave  signal  generators. 

Indeed,  in  1996,  Slonczewski3  and  Berger2  predicted  the 
spin-polarized  phenomenon  in  which  a  polarized  current  can 
exert  a  torque  on  the  magnetization  of  a  ferromagnetic  layer. 
For  large  enough  currents,  this  torque  can  lead  to  switching 
and/or  precession  of  the  magnetization,  see  Fig.  1.  Then  the 
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magnetic  precession  of  the  free  layer  can  lead,  through  the 
GMR  effect,  to  an  oscillating  dipolar  held  in  the  form  of  a 
microwave  voltage  signal  and  turn  the  valve  into  a  spin  tor¬ 
que  nano-oscillator  (STNO).  This  nano-oscillator  is,  in  prin¬ 
ciple,  tunable  over  a  broad  frequency  band,  about  40  GHz,4 


FIG.  1.  Schematic  representation  of  a  spin  torque  nano-oscillator.  This 
nano-oscillator  consists  of  at  least  two  layers  (about  lOOnm  in  lateral  size) 
of  ferromagnetic  materials  separated  by  a  nonmagnetic  material  layer  or 
spacer. 
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which  renders  it  ideal  for  many  applications,  including:  tele¬ 
communications,  e.g.,  wireless  systems;  radar,  e.g.,  air  traffic 
control,  weather  forecasting,  and  navigation  systems.  But  the 
microwave  power  emitted  by  a  single  valve  is  very  small, 
about  1  nW,5  which  is  not  adequate  for  on  chip  applications. 
A  possible  solution  to  this  problem,  which  has  been  proposed 
by  various  groups,6-10  is  to  synchronize  several  STNOs  so 
that  a  coherent  signal  with  a  common  frequency  and  phase 
can  be  extracted  from  the  ensemble  to  produce  a  more 
powerful  (on  the  order  of  micro-watts)  microwave  signal.  To 
date,  there  is  no  report,  however,  that  even  just  five  STNOs 
connected  in  series  can  be  synchronized.  Thus,  an  alternative 
solution  to  achieve  practical  power  is  to  boost  the  power  of  a 
single  STNO.  Indeed,  large  power  (over  1  /(W)  in  single 
STNOs  has  been  recently  demonstrated.1 1  In  this  manuscript, 
we  focus  on  the  first  alternative  as  we  are  motivated  by  all 
previous  works  that  have  tried  to  address  the  problem  of 
synchronization. 

As  noted  by  the  2007  Nobel  Laureate  Professor  Albert 
Fert:  “the  synchronization  of  STNOs  raises  complex  prob¬ 
lems  that  are  new  in  spintronics  and  related  to  the  general 
field  of  Dynamics  of  Nonlinear  Systems.”6  Presumably,  the 
problems  that  Professor  Fert  had  envision  include:  under¬ 
standing  and  classifying  the  various  coherent  states  that  an 
ensemble  of  STNOs  can  produce,  finding  conditions  for  the 
existence  and  stability  of  such  coherent  states,  determining 
the  effects  of  different  couplings  and  connection  topologies, 
establishing  scaling  laws  of  microwave  power  output  for 
large  ID  and  2D  multi-layers,  and  conducting  transformative 
research  to  translate  theory  into  actual  device  realizations  of 
STNOs.  These  problems,  and  many  other  related  issues, 
remain  open  and  are  currently  the  subject  of  intense  research 
efforts,  analytically,  computationally,  and  experimentally, 
across  the  world. 

Thus,  in  2005,  Kaka  and  collaborators  from  the  National 
Institute  of  Standards  and  Technology  (NIST)  reported  in 
Nature,7  the  first  experiments  that  show  that  two  spin  torque 
nano-oscillators  tend  to  phase  lock  when  they  are  in  close 
proximity  of  one  another.  The  coupling  in  this  case  is 
“soft”  as  it  depends  on  the  magnetic  fields  produced  by  each 
nano-oscillator.  Soon  after,  Grollier  et  al.6  investigated,  com¬ 
putationally,  the  behavior  of  a  ID  array  of  N=  10  STNOs 
electrically  coupled  in  series.  Their  study  showed  that  the 
AC  current  produced  by  each  individual  oscillator  can  also 
lead  to  synchronization  and  that,  collectively,  the  microwave 
power  output  of  the  array  is  significantly  larger  than  that  of 
an  individual  valve.  In  a  follow-up  study,  Persson  et  al* 
mapped  out  numerically  the  region  of  synchronization  of  the 
ID  serially  connected  array  considered  by  Grollier  for  the 
special  case  of  N=  2  STNOs.  The  critical  observation,  albeit 
disappointing,  points  out  that  synchronization  is  difficult  to 
achieve  because  the  region  of  parameter  space  where  this  re¬ 
gime  exists  is  rather  small.  Subsequently,  Li  et  al.1'  showed 
that  this  difficulty  was  due,  mainly,  to  the  coexistence  of 
multiple  stable  attractors,  which  suggests  that  the  synchroni¬ 
zation  regime  is  highly  sensitive  to  initial  conditions.  On  a 
more  positive  note,  a  joint  effort  by  researches  from  the 
Army  Research  Laboratory  (ARL)  and  from  NIST  produced 
recently  the  first  demonstration  of  the  ability  of  a  single 


STNO  to  radiate  energy  through  space/  At  about  250  pW 
and  high  frequency  of  9  GHz,  the  generated  signal  carried 
lower  power  than  expected  from  the  previous  theoretical 
studies7  but  it  was  able  to  travel  through  air  over  a  distance 
of  1  m. 

In  the  present  work,  we  provide  a  detailed  description  of 
the  nature  of  the  bifurcations  that  lead  into  and  out  of  the 
synchronization  regime  in  the  ID  array  of  two  serially  con¬ 
nected  STNOs  considered  by  Persson  et  al*  Although 
Persson’s  work,  and  also  our  work,  is  focused  on  only  two 
oscillators,  we  believe  that  a  better  understanding  of  the  na¬ 
ture  of  the  bifurcations  in  this  small  array  can  provide  help¬ 
ful  insight  to  achieve  synchronization  in  larger  arrays.  The 
bifurcation  analysis  shows  bistability  between  in-phase  and 
out-of-phase  patterns  of  oscillations,  which  emerge  via  back- 
to-back  Hopf  bifurcations  from  a  branch  of  nontrivial 
saddle-node  equilibria.  There  are,  in  fact,  two  distinct  pairs 
of  such  limit  cycles.  When  the  applied  magnetic  field  occurs 
in  a  direction  that  leads  to  reflectional  symmetry  in  the  gov¬ 
erning  equations  then  it  might  be  possible  to  manipulate  only 
the  input  current  to  bring  together  the  cycles  until  they  merge 
with  one  another  in  a  gluing  bifurcation.1  ’  But  if  the  applied 
magnetic  field  breaks  the  symmetry  of  the  equations  then 
two  parameters  must  be  varied  to  glue  together  the  cycles. 
This  is  the  case  because  the  gluing  bifurcation  is  a  global 
bifurcation  that  is  facilitated  by  the  presence  of  reflectional 
symmetry.  That  is,  when  reflectional  symmetry  is  present  in 
the  system  the  cycles  can  be  considered  as  mirror  images  of 
one  another.  Thus,  generically,  only  one  parameter  needs  to 
be  varied  to  move  the  two  cycles  and  the  codimension  of  the 
bifurcation  is  one.  In  the  absence  of  reflectional  symmetry, 
each  cycle  can  be  moved  independently  of  one  another  and 
thus  the  bifurcation  becomes  codimension  two.  While  the 
symmetry  of  the  system  is  important,  from  a  mathematical 
standpoint,  for  the  creation  of  the  gluing  bifurcation  the  criti¬ 
cal  observation  that  we  wish  to  emphasize  from  a  physics 
standpoint  is  that  the  bistability  of  the  out-of-phase  oscilla¬ 
tions  is  what  seems  to  make  the  synchronization  state  very 
difficult  to  achieve.  As  the  basins  of  attraction  of  these  two 
solutions  compete  for  stability,  it  appears  that  the  volume  of 
phase  space  that  the  basin  of  attraction  of  the  out-of-phase 
solution  occupies  is  much  larger  than  that  of  the  synchroni¬ 
zation  regime.  The  good  news  is,  however,  that  a  change  in 
the  direction  of  the  applied  field  can  significantly  enlarge  the 
basin  of  attraction  of  the  synchronization  state  and  reverse 
the  preferred  state  towards  synchronization.  This  type  of 
enhanced  synchronization  offers  an  alternative  to  the  use  of 
delays  as  it  was  originally  proposed  by  Persson  et  al 8 
Another  bit  of  good  news  is  that  the  synchronization  regime 
appears  to  be  significantly  robust  to  fluctuations  of  parame¬ 
ters  and,  in  turn,  to  variations  in  the  frequency  of  oscillation 
of  each  individual  spin  valve.  In  particular,  when  the  spin 
valves  are  non-identical,  differences  in  the  anisotropy  and 
demagnetization  fields  and  gyromagnetic  ratio  can  lead  to 
dispersion  among  the  free-evolving  frequency  of  each  indi¬ 
vidual  spin  valve.  However,  careful  tuning  of  the  DC  current 
and  of  the  resistors  in  the  circuit  can  lead  the  array  to  over¬ 
come  the  frequency  mismatch  and  for  the  system  dynamics 
to  evolve  towards  the  synchronization  attractor.  This  result  is 
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expected  because  numerical  computations  show  that  all 
Lyapunov  exponents  transverse  to  the  synchronization  mani¬ 
fold  are  negative;  so  that  small  perturbations  or  fluctuations 
in  system  parameters  can  be  absorbed  by  the  array  while  the 
solution  trajectories  asymptotically  converge  to  the  synchro¬ 
nization  state.  From  the  standpoint  of  the  bifurcation  analy¬ 
sis,  we  also  wish  to  emphasize  that  differences  in  parameter 
values  tend  to  shift  the  onset  of  oscillations  shown  in  the 
computational  bifurcation  diagrams  while  the  overall  struc¬ 
ture  of  the  diagrams  is  preserved.  We  do  not  attempt  in  this 
work  to  quantify  the  frequency  mismatch  nor  the  shift  in  the 
location  of  the  Hopf  bifurcations  that  lead  to  limit  cycle 
oscillations.  Those  tasks,  and  a  more  comprehensive  analysis 
of  the  effects  of  the  nonhomogeneity  of  parameter  values 
and  an  investigation  of  stochastic  effects,  such  as  noise  in 
the  system,  are  part  of  ongoing  work.  We  expect  to  report 
those  results  in  a  follow-up  manuscript. 

II.  MODELING 

A.  Single  STNO  dynamics 

An  originally  unpolarized  electric  current  /,  in  units  of 
Amp,  is  applied  to  a  fixed  magnetic  layer  whose  magnetiza¬ 
tion  is  represented  by  M.  As  the  electrons  pass  through  the 
layer,  their  spins  move  (and  flip  if  necessary)  to  align  their 
orientations  to  that  of  the  fixed  layer,  thus,  creating  a  spin- 
polarized  current,  see  Fig.  1  for  a  schematic  diagram. 

The  electrical  potential  that  exists  across  the  nonmag¬ 
netic  layer  (labeled  spacer)  allows  the  spin-polarized  current 
to  preserve  its  polarization.  So  when  the  spin-polarized  cur¬ 
rent  enters  the  free  magnetic  layer,  it  exerts  a  torque  on  its 
magnetization  in.  According  to  Newton’s  third  law,  the 
amount  of  torque  is  directly  proportional  (and  of  opposite 
sign)  to  the  difference  in  the  magnetization  of  the  spins  in 
the  polarized  current  and  those  of  the  free  layer.  We  will 
assume  the  layers  to  be  uniform  so  that  the  spin  precession  is 
proportional  to  —in  x  / /err,  where  He ff  is  the  effective  mag¬ 
netic  field,  which  consists  of  the  exchange  field,  //exchange, 
among  magnetic  moments,  a  surface  anisotropy  field, 
H anisotropy,  which  defines  a  preferred  direction  of  magnetiza¬ 
tion,  a  demagnetization  field  H demagnetization,  and  the  applied 
magnetic  field  H  applied-  Collectively,  the  effective  field 
becomes 


where  y  is  the  gyromagnetic  ratio,  in  units  of  „ '  ,  where  ns 
represents  nanoseconds,  while  X  serves  as  the  magnitude  of 
the  damping  term,  in  dimensionless  units.  In  the  spin  torque 
term,  a  has  a  unit  of  Oe  and  is  proportional  to  the  electrical 
current  density14,18'19  which  can  be  written  as  a  = 
where  S0  is  the  constant  magnitude  of  the  average  magnet¬ 
ization  vector  S(t).  in  units  of  Oe,  so  that  in  =  S/So  is  the 
dimensionless  unit  vector  in  the  direction  of  .S',  g  is  a  function 
of  the  polarization  factor  P,  which  will  be  assumed  to  be 
exactly  one  in  dimensionless  units,  h  =  6.582  x  1CT16  is 
Planck’s  constant  in  units  of  eV  s,  V  =  3.0732  is  volume  in 
units  of  cm3,  e  =  1.602  x  1 0  1 9  is  the  elementary  charge  in 
units  of  Coulombs. 

B.  Series  array  of  STNOs 

We  now  consider  a  ID  array  of  STNOs  connected  in  se¬ 
ries,  as  is  shown  in  Fig.  2,  and  derive  the  governing  equa¬ 
tions  for  the  general  case  of  N  STNOs  though  the  remaining 
work  is  focused  on  the  particular  case  of  A  =  2  studied  by 
Persson  et  al ,8 

Following  the  work  of  Grollier  et  al.,6  we  assume  the 
standard  equation  for  the  resistance  (in  units  of  Q)  of  the  zth 
oscillator  to  be  /?,(/)  =  Rqi  —  A/?,-cos0,-(f),  where  0,(t)  is  the 
angle  between  the  magnetization  of  the  fixed  and  free  ferro¬ 
magnetic  layers,  Rq,  is  the  mean  while  A/?,  is  half  the  differ¬ 
ence  between  the  resistances  in  the  parallel,  RPi,  and  the 
anti-parallel,  Rapi,  magnetization  states,  respectively.  That 
is,  R ( ),  =  ( Rapi  +  Rp,)/2  and  A Rj  =  (Rapi  —  R  in) /2.  The 
input  /0  is  a  known  DC  current.  To  determine  the  instantane¬ 
ous  current  through  the  yth  STNO  element,  we  combine 
Kirchoff’s  Current  Law  and  Ohm’s  Law  to  produce  a  simple 
current  divider  equation: 


J2R‘+Rc 

;=  i 

Because  the  right-hand  side  of  Eq.  (2)  is  independent  of 
j,  the  current  must  be  the  same  in  all  oscillators.  Removing 
the  j  index  and  substituting  R,  into  Eq.  (2)  produces,  after 
some  manipulation,  the  following  equation  for  the  current: 


Z/eff  —  //exchange  T  //anisotropy  T  //demagnetization  T  //applied- 

Also  it  can  be  shown  that  the  spin-transfer  torque  is  pro¬ 
portional  to  in  x  (in  x  M).  Energy  dissipation  effects  such  as 
those  due  to  spin  scattering  lead  to  a  damping  term  proportional 
to  in  x  '-/j/.  Together,  these  quantities  govern  the  dynamics  of 
the  free  magnetization  layer  through  the  Landau-Lifshitz- 
Gilbert-Slonczewski  (LLGS)2’14-17  equation: 

damping 

precession  ^ 

dm  '  *  „  dirt 

—  =  —  y  m  x  H,  ff  +  a  m  x  — 
dt  dt 

spin  transfer  torque 

- - : - 

—  yag(P,m  ■  M)ih  x  (m  x  M),  (1) 


FIG.  2.  Circuit  array  of  spin  torque  nano-oscillators  connected  in  series. 
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N 

k<  -  'y  R0i 

m  ^ - /o.  o) 

N  [  A Rj  cos  6j(t) 


Rc+y:  Roi 

i=  1 


Notice  that  the  numerator  is  time  invariant.  In  fact,  this 
numerator  is  a  good  approximation  for  the  DC  current  in  the 
oscillator  circuit  branch.  Expanding  Eq.  (3)  in  a  first  order 
Taylor  approximation,  we  can  rewrite  I(t)  in  the  simplified 
form: 


where 

l  DC  = 


I(t)  —  IdC  fijCOsfl^f)^  ) 


Rc  ,  A  0  ARi 

— - 10  and  0Mi  = - — 


(4) 


Rc  +  ^>2  Roi 


Rc  +  ^ 2  R°‘ 


We  note  that  the  current  /  appears  in  the  spin  torque 
term  of  the  LLGS  Eq.  (1)  through  the  parameter  a.  Thus, 
assuming  a  polarization  factor  g=  1,  we  arrive  at  the  follow¬ 
ing  model  for  the  array  of  N  STNOs  electrically  coupled  in 
series 


dm j 
dt 


-ymj  x  Heff  +  hhj  x 


dmj 

dt 


-y 


2  S0Ve 


(  \ 

I  DC  (H-  ^ar/cos0/(OJ«C  x  x  44). 

1=1 


(5) 


It  is  observed  from  Eq.  (5)  that  the  series  configuration 
results  in  an  all-to-all  coupling  scheme  where  each  individ¬ 
ual  STNO  is  influenced  by  the  angles  between  the  free  and 
fixed  layers  of  every  other  STNO.  Thus,  the  symmetries  of 


the  series  array  of  N  STNOs  are  described  by  the  group  Sn, 
which  is  the  group  of  all  permutations  of  N  objects.  While 
these  equations  are  valid  for  up  to  j  =  I...N  oscillators,  our 
aim  in  this  work  is  to  focus  on  the  case  N  =  2  considered  by 
Persson  et  al.8  But  first  we  conduct  computer  simulations  of 
the  governing  equations  (5)  in  order  to  get  insight  into 
the  type  of  transitions  that  lead  to  the  synchronization  state. 
We  employ  the  relations  introduced  by  Murugesh  and 
Lakshmanan1814  for  the  following  terms:  //anisotropy 
=  K(m  ■  e\\)e\\*  where  k  =  45  Oe,  ey  =  [sin0||Cos</>||,  sinfly 

sin</>y,COS0y]r  is  dimensionless,  //demagnetization  =  -4tt50 
{N\m\X  +  A?2»J2,V  +  where  Nh  i=  1,2,3  are 

dimensionless  constants  with  N\  +  N2  +  IV3  =  1,  and 
{x,y,  z}  are  the  orthonormal  unit  vectors.  However,  we  devi¬ 
ate  slightly  in  considering  the  applied  magnetic  field  to  lie 
on  the  vz-plane  instead  of  the  z-axis,  so  that  //applied 
=  ha  [0,  sin0/,,  cos0/,]r,  where  ha  is  in  units  of  Oe.  In  what 
follows  we  assume  the  direction  of  demagnetization  to  be 
along  the  i-axis  so  that  N  =  j  1 . 0, 0|7  ,  thus,  creating  a  yz- 
plane.  We  also  assume  fly  =  0  so  that  ey  =  [0,0, 1],  which 
produces  an  easy  axis  in  the  z-direction.  Finally,  we  assume 
the  direction  of  polarization  of  the  spin-polarized  current  to 
remain  constant  along  the  z-direction  M  =  z. 

In  the  computational  work  of  Persson  et  al.,8  it  was 
reported  that  the  magnetization  states  m.j{t)  relax  to  stable 
equilibrium  states  for  small  and  for  very  large  positive  values 
of  the  input  DC  current,  IDC.  It  was  also  reported  that  non- 
synchronized  oscillations  occur  for  most  intermediate  values 
of  IDC  while  synchronized  oscillations  are  rare — as  they  occur 
on  very  small  regions  neighboring  the  equilibrium  states. 
Something  very  similar  occurs  in  our  simulations,  see  Figure 
3,  except  that  now  the  large  values  of  current  where  equilib¬ 
rium  states  appear  are  negative.  This  inversion  of  sign  is  due 
to  the  fact  that  in  our  formulation  of  the  LLGS  Eq.  (1)  we 
have  employed  the  notation  introduced  by  Lakshmanan,1' 
which  contains  a  negative  sign  in  front  of  the  gyromagnetic 
ratio  y  as  oppose  to  the  positive  sign  used  by  Persson  et  al.  Up 
to  a  conjugacy  in  sign,  these  two  formulations  are  equivalent 
So  we  can  proceed  with  the  exploration  of  the  dynamics.  As 
reported  by  Persson  et  al.,  our  simulations  show  that  for 


FIG.  3.  Collective  behavior  of  two  STNOs  coupled  in  a  series  array  through  an  external  electrical  current  IDC.  Parameters  are:  y  =  0.0 1 76  0)rls-  A  =  0.008, 
011  =  0  rad,  0h  =  0  rad,  ha  =  300  Oe,  K  =  45  Oe,  S0  =  8400/(4ji)  Oe,  N  =  [1 , 0,  Of ,  Roi  =  0. 1  Q,  Rc  =  50  Q,  and  A R,  =  0.03  Q. 
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intermediate  values  of  IDC  the  dynamics  is  attracted  to  stable 
limit  cycle  oscillations.  There  are,  indeed,  two  types  of  limit 
cycles.  One  of  them  corresponds  to  the  out-of-phase  oscilla¬ 
tions,  OP,  where  the  two  magnetization  states  oscillate  with 
the  same  waveform  and  amplitude  but  out  of  phase  by  half  a 
period.  The  other  cycle  corresponds  to  the  in-phase  oscilla¬ 
tions,  IP,  leading  to  a  complete  synchronization  state.  Next  we 
study  in  more  detail  the  nature  of  the  bifurcations  that  lead 
into/out-of  the  synchronization  state.  We  would  like  to 
emphasize  that  such  bifurcation  study  was  not  discussed  in 
Persson  et  al 8  but  it  can  be  an  important  tool  to  help  us  look 
for  clues  to  increase  the  region  of  stability  of  the  synchronized 
pattern  of  collective  behavior  with  the  ultimate  goal  of  achiev¬ 
ing  higher  power  output  through  an  array  of  spin  valves. 

III.  COMPUTATIONAL  BIFURCATION  ANALYSIS 

In  order  to  understand  the  nature  of  the  bifurcations  that 
lead  into  and  out  of  the  synchronization  state,  we  perform 
next  a  computational  bifurcation  analysis,  with  the  aid  of  the 
continuation  software  package  auto,20  of  the  governing 
equations.  For  convenience,  we  convert  to  complex  stereo¬ 
graphic  coordinates18'14  through  the  following  change  of 
variables: 


iriji  +  imJ2 
1  +  mj3 


(6) 


Direct  and  tedious  calculations  using  Eq.  (6)  lead  to  the 
following  set  of  differential  equations  in  stereographic 
coordinates: 


(1 


i-7  yha2  / ,  2\ 

-  yazj  +  lyhaizj  +  -y-  ( 1  +  zj) 


+  im\\K.y 


cos0||Z2  —  —  sin0|  |  (e" 


iyAnS0 

(i  +  M2) 


N3(l  -  \zj\2)zj 


Nj_ 

2 

(N 


(1 


-z2 

J 

N2) 
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As  a  first  case,  we  set  0/,  =  0  so  that  the  applied  mag¬ 
netic  field  //applied  lies  exactly  along  the  z-axis  as  it  was  origi¬ 
nally  considered  by  Persson  et  al.8  We  then  use  Eq.  (7)  to 
generate  the  one-parameter  bifurcation  diagram  in  Fig.  4, 
which  illustrates  how  the  magnetization  state  changes  as  a 
function  of  input  current.  As  expected,  for  positive,  small 
and  large  values  of  the  input  current,  the  magnetization  state 
relaxes  to  the  stable  trivial  or  zero  equilibrium  state,  which  is 
indicated  as  a  solid  (red)  line.  For  large  negative  values  of 
the  input  current,  in  particular,  two  branches  of  nontrivial 
stable  equilibrium  points  emerge  via  the  saddle  node  bifurca¬ 
tions,  labeled  SN\  and  ,S'/V2.  Along  these  branches,  the  spin 
valves  cannot  oscillate,  as  they  remain  in  a  steady  state. 
Many  other  equilibrium  points  also  exist  over  a  wide  range 
of  values  of  IDC  but  they  are  mostly  unstable  (dashed  line). 
As  the  input  current  increases  from  the  saddle-node 


bifurcation  points  SNX,  the  nontrivial  equilibrium  point 
remains  almost  constant  until  it  loses  stability  and  two 
branches  of  periodic  solutions  then  emerge  via  back-to-back 
Hopf  bifurcations,  labeled  HB4  and  HB5.  One  branch  (blue) 
corresponds  to  the  out-of-phase  pattern,  OP,  and  the  other 
one  (green)  corresponds  to  the  complete  synchronization 
state,  IP.  It  is  not  surprise  that  these  are  the  only  oscillations 
that  appear  for  they  are  the  most  common  patterns  of  collec¬ 
tive  behavior  observed  in  a  system  of  two  identical  oscilla¬ 
tors  coupled  symmetrically,  as  it  the  case  of  our  STNO 
array.  Furthermore,  it  can  be  shown  that  the  oscillations 
emerge  through  symmetry-breaking  bifurcations  in  which 
the  synchronous  state  preserves  the  S2  symmetry  of  the  array 
while  the  out-of-phase  pattern  breaks  it.  Now,  since  the  mag¬ 
netic  field  was  applied  directly  along  the  z-axis,  the  govern¬ 
ing  equations  exhibit  also  a  reflectional  Z2  symmetry  for  the 
internal  dynamics  of  each  individual  spin  valve,  which  mani¬ 
fests  as  a  mirror-image  symmetry  across  the  X\  =  0  axis  in 
the  bifurcation  diagram.  This  symmetry  leads  to  a  second 
pair  of  limit  cycle  oscillations  that  emerge  from  the  back-to- 
back  Hopf  points  HB2  and  //fi3  with  the  same  characteristics 
as  those  of  HB4  and  HB5.  Now  as  the  coupling  DC  current 
IDC  increases  the  amplitude  of  the  limit  cycles  in  each  branch 
increases  gradually  until  the  cycles  merge  with  one  another 
in  a  gluing  bifurcation  near  Idc  =  0.  The  onset  of  the  gluing 
bifurcation  occurs  just  before  the  Hopf  bifurcation  HBX, 
which  appears  off  of  the  trivial  equilibrium. 

Additionally,  we  wish  to  emphasize  that  when  the  spin 
valves  are  non-identical,  fluctuations  in  system  parameters, 
such  as  anisotropy  and  demagnetization  fields  and  gyromag- 
netic  ratio,  can  shift  the  onset  of  the  Hopf  bifurcation  points 
that  lead  to  limit  cycle  oscillations.  However,  the  overall 
structure  of  the  bifurcation  diagrams  remains  the  same  pro¬ 
vided  that  the  fluctuations  are  not  extremely  large.  In  a 
follow-up  manuscript,  we  plan  to  explore  the  effects  of  the 
nonhomogeneity  of  parameter  values  as  well  as  the  effects  of 
noise  on  the  bifurcation  structure. 

Up  to  now,  we  have  described  somehow  mathematically 
complex  behaviors  associated  with  the  bifurcation  diagrams, 
as  illustrated  in  Fig.  4,  where  IP  and  OP  oscillations  overlap 
across  a  wide  range  of  IDC  values.  From  an  experimental 
physics  standpoint,  our  main  interest  is  to  better  understand 
the  conditions  for  the  existence  and  stability  properties  of  the 
synchronization  behavior,  as  it  can  provide  benefits  in  oper¬ 
ating  an  array  of  STNOs  to  increase  the  overall  radiative 
power.  The  regions  of  existence  are  outlined  through  the 
bifurcation  diagrams  while  the  stability  properties  can  be 
studied  through  many  different  approaches.  One  such 
approach  is  to  compute  the  transverse  Lyapunov  exponents 
(TLE)  of  the  synchronization  manifold-1'"2  of  the  coupled 
system.  To  do  this,  we  start  by  re-scaling  time  so  that  t 
=  1  and  replace  the  complex  stereographic  coordinates 

zj,  j  =  1 , 2,  with  real  coordinates  (xj.  yj),  where  zj  =  Xj  +  yji, 
so  that  the  governing  Eq.  (7)  becomes 

=fi.Xj,yj)  +  hl(xuyux2,y2),  g 

jj  =  g(xj,yj)  +  h2(xi,yi,X2,y2), 


where 
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FIG.  4.  (Top)  One  parameter  bifurcation 
diagram  of  changes  in  magnetization  in 
a  ID  array  of  N  —  2  STNOs  connected  in 
series.  Blue  circles  indicate  out-of-phase 
oscillations  while  green  circles  indicate 
synchronized  limit  cycle  oscillations. 
Filled-in  (empty)  circles  indicate  stable 
(unstable)  oscillations.  (Middle) 
Visualization  of  the  dynamics  on  the  unit 
sphere  and  (bottom)  on  the  stereographic 
plane.  Parameters  are:  7  =  0.0176q^, 
A  =  0.008,  0||  —  0  rad,  0/7  =  0  rad,  ha 
=  300  Oe,  k  =  45  Oe,  S0  8400/(4tt) 
Oe,  N  =  [1,0, 0]r,  R0i  =  0.10,  Rc  =  50 
Q,  and  AS,  =  0.03  Q. 
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represent  the  internal  dynamics  of  the  Xj  and  yj  components 
of  the  /th  nano-oscillator,  respectively,  and  rj  =  xj  +  yj.  The 


terms  hj  are  the  (all-to-all)  coupling  functions  that  connect 
the  two  nano-oscillators  together,  they  are  given  by 

hi  (xj,yj)  =  -axj+aAyj, 
h2(xj,yj )  =  -afocj  -  ay}. 


The  coupling  strength  is  set  by  the  parameter  a.  If  we  assume 
that  all  the  oscillators  have  the  same  resistance  =  Rq  and 
magnetoresistance  ARj  =  A R  then  the  coupling  strength  for 
N  =  2  can  be  written  as 


a  = 


Ti 

2SQVelDC 
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Let  Uj  =  ( Xj,yj),  j  =  1,2.  We  now  transform  Eq.  (8)  to 
the  transversal  coordinates  x±  =  u\  —  u2  and  consider  only 
small  perturbations  to  the  synchronization  manifold,  i.e., 
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Mi  «  ui  .  Then  the  dynamics  along  the  transverse  direction  to 
the  synchronization  manifold  becomes 


*l,  =  /0 i , yi)  -fix2,yi)  -  ax±1  +  aXx±2 , 

Tj-2  =  gixuyi)  -  g(x 2,yi)  -  ahiU  -  ax±1. 

We  can  now  proceed  in  two  ways.  Expand /(X2.V2)  in  a 
Taylor’s  series  expansion  about  (xi,yi)  or  expand  f(x\ , yi ) 
about  the  second  coordinate.  We  choose  the  former  approach 
for  both  f  and  g,  and  also  for  the  terms  in  a.  Neglecting 
higher  order  terms  of  x± ,  we  obtain,  in  matrix  form,  the  lin¬ 
earization  of  Eq.  (8)  transverse  to  the  synchronization 
manifold 


x±  =  (J  +  K)x±, 


(10) 


where  J  is  the  Jacobian  matrix  of  the  vector  field  F 
=  (f(x,y),g(x,y))  evaluated  at  the  synchronization  manifold, 
i.e.,  J  =  dF(Xl0ll),  where  xs  =  xi  =,  ...,xN,  ys  =  yi  =,  ...,yN, 
and  K  is  the  matrix  that  results  from  the  linearization  of  the 
coupling  terms 
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The  synchronization  state  u\  =  U2  is  said  to  be  asymp¬ 
totically  stable  if  x±  — »  0  as  t  — *  00  or  if  the  transverse 
Lyapunov  exponents,  which  are  the  Lyapunov  exponents 
associated  with  the  linearized  equations  (10),  are  all  nega¬ 
tive.  Numerical  computations,  see  Fig.  5,  reveal  that  the  sum 
of  TLEs  is  always  negative  except  for  the  interval 
1  DC  £  [-1700,-1 660] ,  which  is  enhanced  in  the  inset. 
Recall  from  the  one-parameter  bifurcation  diagram  of  Fig.  4 
that  the  subinterval  Iqc  £  [—1700,  —1675]  is  the  region 
between  the  Hopf  points  HB4  and  IIB5  where  only  the  out- 
of-phase  solution  exists,  so  the  calculation  of  TLEs  picks  up 
the  instability  of  the  unstable  non-trivial  steady  state  and  that 


FIG.  5.  Transverse  Lyapunov  exponents  computed  as  a  function  of  input 
current  for  an  array  of  N  —  2  STNOs  connected  in  series  with  applied  field 
along  the  z-direction,  i.e..  Oh  =  0  rad.  All  other  parameters  are  the  same  as 
in  Fig.  4. 


explains  why  the  sum  of  TLEs  is  positive.  As  IDC  reaches 
the  right-end  of  that  subinterval,  the  synchronized  solution 
emerges  via  a  subcritical  Hopf  bifurcation  at  HB 5,  so  the 
sum  of  TLEs  is  still  negative.  However,  as  IDC  increases  fur¬ 
ther  within  the  subinterval  Iqc  £  [—1675,  —1660],  the 
synchronized  pattern  becomes  less  unstable  and  the  sum  of 
TLEs  decreases  until  it  eventually  crosses  zero  near 
I dc  =  —1660,  at  which  point  the  synchronized  state  becomes 
stable.  From  then  on,  the  sum  of  TLEs  decreases  monotoni- 
cally  while  the  synchronized  pattern  increases  stability  until 
the  gluing  bifurcation  point  is  reached  near  Iqc  ~  —100. 
Past  this  point,  the  sum  of  TLEs  increases  until  IDC  reaches 
the  Hopf  point  HBX.  To  the  right  of  HBX  the  trivial  equilib¬ 
rium  is  stable  and  so  the  sum  of  TLEs  is  negative.  Similarly, 
in  the  region  where  the  sum  of  TLEs  is  negative,  all  individ¬ 
ual  TLEs  are  also  negative;  though,  they  are  not  shown  for 
brevity.  Consequently,  in  that  region  there  is  a  contraction  of 
the  phase  space  volume  of  the  dynamics  that  is  transversal  to 
the  synchronization  manifold.23  Negative  TLEs  also  account 
for  the  robustness  exhibited  by  the  synchronization  state  to 
fluctuations  in  parameters.  Indeed,  we  have  verified  through 
computer  simulations  that  when  the  spin  valves  are  non¬ 
identical,  changes  in  the  anisotropy  and  demagnetization 
fields  and  gyromagnetic  ratio  can  lead  to  dispersion  among 
the  free-evolving  frequency  of  each  individual  spin  valve. 
Small  frequency  dispersions  can  still  lead  the  system  dynam¬ 
ics  towards  the  synchronization  attractor  by  careful  tuning  of 
the  DC  current  and/or  the  resistors  in  the  circuit  which  in 
turn  control  the  parameter  AGMR.  We  do  not  attempt,  how¬ 
ever,  to  measure  the  largest  dispersion  that  can  be  overcome 
by  the  system.  That  task  is  part  of  ongoing  work. 

Now,  since  the  nano-oscillators  are  mutually  coupled  in 
an  all-to-all  fashion,  the  synchronized  state  may  not  neces¬ 
sarily  be  globally  asymptotically  stable  and  may  still  depend 
on  initial  conditions.  Indeed,  we  already  know  from  the 
bifurcation  diagram  that  the  out-of-phase  solution  coexists 
and  is  stable  within  the  same  parameter  region  of  the  syn¬ 
chronization  state.  Each  one  of  these  solutions  has  its  own 
basin  of  attraction.  We  do  not  attempt  in  this  work  to  visual¬ 
ize  those  basins  of  attraction.  However,  we  try  to  get  a  sense 
of  their  size  by  measuring  the  coherence  parameter,  also 
known  as  the  Kuramoto24'23  order  parameter 


r 


1 

N 


N 


k=  l 


0  <  r  <  1. 


In  Fig.  6,  we  compute  the  Kuramoto24,25  order  parameter 
r  when  Of,  =  0  over  a  region  where  one  of  the  pairs  of  limit 
cycles  coexist  as  a  function  of  DC  current.  This  parameter 
serves  as  a  measure  of  coherence  of  the  phase  dynamics  of  an 
ensemble  of  oscillators.  Indeed,  direct  calculations  show  that 
when  ;•  =  0  the  phase  dynamics  is  asynchronous  while  r  =  1 
corresponds  to  full  or  complete  synchronization.  Intermediate 
values  of  r  can  be  associated  with  more  complicated  states 
where  the  oscillators  organize  themselves  into  clusters  of  in- 
phase,  traveling  waves,  and  other  related  patterns.  Thus,  the 
calculation  of  the  coherence  parameter  suggests  that  when 
dh  =  0  the  coupled  system  is  highly  sensitive  to  changes  in 
input  current.  Sometimes  the  phase  dynamics  of  the  two  nano- 
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FIG.  6.  Coherence  parameter  r  computed  as  a  function  of  input  current  for  an 
array  of  N  =  2  STNOs  connected  in  series  with  applied  field  along  the  z-direc- 
tion,  i.e.,  0/,  =  0  rad.  r  =  1  is  associated  with  full  synchronization  while  r  —  0 
is  indicative  of  asynchronous  behavior.  Parameters  are  the  same  as  in  Fig.  4. 


oscillators  synchronize  but  then  a  slight  change  in  the  input 
current  can  easily  push  the  system  in  and  out  of  the  synchroni¬ 
zation  state.  The  transitions  from  one  pattern  of  oscillation  to 
another  do  not  occur,  in  general,  instantaneously  so  r  can  fluc¬ 
tuate  slightly  between  zero  and  one.  Presumably  we  could  cor¬ 
rect  those  fluctuations  by  integrating  for  longer  periods  of 
time.  However,  we  would  like  to  point  out  that  we  are  already 
integrating  the  equations  for  a  vey  long  time  interval,  so  fur¬ 
ther  increases  can  yield  minor  changes  in  r.  More  importantly, 
the  lack  of  large  open  intervals  where  r  could  remain  close  to 
one,  see  Fig.  6,  suggest  that  the  basin  of  attraction  of  the 
synchronized  oscillations  occupies  a  relatively  small  volume 
of  phase  space.  This,  again,  confinns  the  observation  by 
Persson  et  al.':  that  “nonsynchronized  precession  largely  dom¬ 
inates  the  phase  diagram.”  Similar  results  are  obtained  for  the 
other  pair  of  limit  cycles. 

Figure  7  shows  the  two-parameter  continuation  of  the 
Hopf  bifurcation  points  HB\  —  HB5  and  saddle-node  points, 
in  parameter  space  k  vs  Idc.  Recall  that  k  is  a  coefficient 
associated  with  the  anisotropy  field.  For  small  values  of  k, 
the  Hopf  points  HB2,  ...,HB^  are  significantly  separated 
from  the  saddle-node  points  SN\  and  .S'/V2,  as  is  the  case  in 
the  one-parameter  bifurcation  diagram  of  Fig.  4.  As  k 
increases,  however,  the  locus  of  the  back-to-back  Hopf 
points  move  closer  to  that  of  the  saddle-nodes,  which  remain 
almost  constant.  At  the  same  time,  the  region  of  stable 
synchronized  limit  cycle  solutions,  shaded  gray,  starts  to 
shrink.  For  significantly  larger  values  of  k,  the  locus  of  the 
Hopf  point  has  merged  with  the  locus  of  the  saddle-node 
points  and  the  region  of  stable  limit  cycle  in-phase  oscilla¬ 
tions  has  all  but  disappeared.  For  completeness  purposes,  the 
bifurcation  diagram  is  carried  out  over  negative  values  of  k. 
Observe  that  the  loci  of  the  Hopf  points  HB 2, . ...  IIB4  meets 
the  almost  circular  curve  that  defines  the  loci  of  pitchfork 
bifurcations  of  nontrivial  equilibrium  points.  These  equilib¬ 
rium  states  are  the  continuation  of  the  saddle-node  points  for 
small  negative  values  of  k,  so  they  do  not  show  up  in  Fig.  4. 


FIG.  7.  Two  parameter  bifurcation  diagram  depicts  the  boundary  curves,  as 
functions  of  input  current  and  anisotropy  coefficient,  which  separate  equilibrium 
states  from  limit  cycle  oscillations  in  a  ID  array  of  N=2  STNOs  connected  in 
series.  IP/OP  indicate  in-plane  and  out-of-plane  oscillations,  respectively. 
Parameters  are  the  same  as  in  Fig.  4  but  we  emphasize  0/,  —  0  rad. 

Applying  the  magnetic  field  with  a  small  angle  along  the 
y-di  recti  on,  e.g.,  ()/,  =  n/ 4,  breaks  the  Z2-symmetry  of  the 
governing  equations  but  the  one-  and  two-parameter  bifurca¬ 
tion  diagrams  exhibit  similar  characteristics  to  the  perfectly 
symmetric  case,  i.e.,  0/,  =  0,  as  is  shown  in  Fig.  8(left).  That 
is,  there  are  two  pairs  of  branches  of  periodic  oscillations, 
each  pair  contains  a  branch  of  synchronized  limit  cycle  oscil¬ 
lations  and  a  branch  of  out-of-phase  oscillations.  The  two 
pairs  might  merge  together  again  in  a  gluing  bifurcation. 
However,  in  the  absence  of  the  Z2-symmetry,  it  may  be  nec¬ 
essary  to  vary  two  parameters  simultaneously  for  the  gluing 
of  the  limit  cycles  to  occur.  This  result  shows  that  the  reflec- 
tional  Z2  symmetry  is  important  to  generate  the  gluing  bifur¬ 
cation  via  a  codimension  one  bifurcation;  otherwise  the 
codimension  is  two.  As  we  discussed  earlier,  the  two  patterns 
of  oscillations  emerge  again  via  symmetry-breaking  bifurca¬ 
tions  of  the  network,  and  in  which  the  synchronized  state  pre¬ 
serves  the  S2  symmetry  of  the  array  while  the  out-of-phase 
state  breaks  it.  Now  a  critical  observation  is  that  if  the  angle 
of  the  applied  field  is  further  increased,  e.g.,  0 /,  =  3n/4,  then 
it  is  possible  for  the  IP  and  OP  branches  of  oscillations  to 
switch  position,  as  is  shown  in  the  one -parameter  bifurcation 
of  Fig.  8(right).  This  is  a  critical  observation  because  the 
switch  leads  to  an  interval  of  the  DC  current  where  only  the 
synchronized  state  is  stable  while  the  OP  solution  is  unstable. 
Thus,  in  principle,  on  this  interval  it  should  be  easier  to 
achieve  synchronization.  Furthermore,  we  can  see  in  Fig. 
9(left)  that,  when  0/,  =  n/4,  small  open  intervals  where  r  w 
1  start  to  appear  throughout  the  region  where  the  two  limit 
cycles,  IP  and  OP,  coexist.  This  suggests  that  a  slight  change 
in  the  direction  of  the  applied  field  can  facilitate  the  existence 
of  the  synchronization  state.  Indeed,  increasing  0/,  further  to 
37t/4  shows,  see  Fig.  9(right),  a  very  large  open  interval  of 
DC  current  where  the  coherence  parameter  r  is,  approxi¬ 
mately,  equal  to  one.  Thus,  changing  the  direction  of  the 
applied  field  appears  to  have  the  net  effect  of  changing  the 
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(a)One  parameter  bifurcation  diagram  with  Oh  =  f 


(b)One  parameter  bifurcation  diagram  with  Oh  = 


x  io4 


(c)Two  parameter  bifurcation  diagram  with  Qh  =  \  (d)Two  parameter  bifurcation  diagram  with  Oh  =  x 


FIG.  8.  One-  and  two-parameter  bifurcation  diagrams  indicating  regions  of  existence  and  stability  of  equilibrium  states  and  limit  cycle  oscillations  for  a  ID 
array  of'  /V  —  2  STNOs  connected  in  series,  with  0/,  =  n/A  rad  and  0/,  =  3 tt /  4  rad.  All  other  parameters  are  the  same  as  in  Fig.  4. 


FIG.  9.  Coherence  parameter  r  computed  as  a  function  of  input  current  for  an  array  of  N  —  2  STNOs  comiected  in  series  with  applied  field  defined  by  (left) 
0/,  =  n/A  rad  and  (right)  0/,  —  3tt/4  rad.  All  other  parameters  are  the  same  as  in  Fig.  4. 


basin  of  attractions  so  that  the  one  associated  with  the  IP  so¬ 
lution  now  dominates  a  larger  volume  of  phase  space.  This 
change  is  already  a  significant  improvement  in  achieving 
synchronization  especially  when  we  compare  Figs.  6  and  9. 

IV.  CONCLUSION 

In  this  work,  we  have  studied  the  collective  behavior  of 
an  array  of  spin  torque  nano-oscillators  connected  in  series. 


The  derivation  of  the  governing  equations  for  the  array  leads 
to  an  all-to-all  SN  symmetric  coupled  system  of  ordinary  dif¬ 
ferential  equations,  in  which  each  individual  nano-oscillator  is 
modeled  through  the  Landau-Lifshitz-Gilbert-Slonczewski 
equation.  A  computational  bifurcation  analysis,  focused  on 
two  nano-oscillators,  shows  transitions  between  equilibrium 
states  and  limit  cycles  oscillations  in  a  manner  that  is  consist¬ 
ent  with  related  studies  conducted  by  other  researchers.  In  par¬ 
ticular,  the  bifurcation  diagrams  reveal  two  patterns  of 
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collective  behavior:  an  in-phase  pattern,  in  which  all  oscilla¬ 
tors  are  fully  synchronized  with  one  another,  and  an  out-of- 
phase  pattern.  These  two  patterns  are  generic,  i.e.,  most 
common,  solutions  in  arrays  of  two  identical  oscillators 
coupled  symmetrically.  Furthermore,  the  presence  of  reflec- 
tional  symmetry  in  the  internal  dynamics  of  each  nano¬ 
oscillator  leads  naturally  to  gluing  bifurcations  of  these  two 
cycles.  However,  of  greater  interest  from  a  physics  standpoint 
is  the  region  of  existence  and  stability  properties  of  the  in- 
phase  pattern  because  synchronized  oscillations  can,  in  princi¬ 
ple,  yield  higher  microwave  power  output.  Previous  works 
have  shown  that  such  pattern  is  rather  difficult  to  realize.  The 
work  conducted  in  this  manuscript  provides  answers  to  cir¬ 
cumvent  this  difficulty.  More  specifically,  the  bifurcation 
analysis  and  measures  of  coherence  presented  in  this  manu¬ 
script  reveal  that  a  change  in  the  direction  of  the  applied  held 
can  significantly  enhance  the  basin  of  attraction  of  the  syn¬ 
chronous  state  without  the  need  of  adding  delay  into  the 
model  equations.  Work  in  progress  includes  visualization  of 
the  basins  of  attraction  and  a  similar  analysis  of  the  collective 
behavior  of  larger  arrays,  which  we  can  expect  to  show  more 
complicated  patterns  of  oscillations.  Additionally,  we  plan  to 
carry  out  a  comprehensive  analysis  of  the  robustness  of  the 
system  to  fluctuations  in  parameters,  quantification  of  the  non¬ 
homogeneity  of  parameter  values,  and  an  investigation  of  sto¬ 
chastic  effects  such  as  noise  in  the  system.  We  expect  to 
report  those  results  in  a  follow-up  manuscript. 
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